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Abstract 

Landscape features of anthropogenic or natural origin can influence organisms' 
dispersal patterns and the connectivity of populations. Understanding these 
relationships is of broad interest in ecology and evolutionary biology and pro- 
vides key insights for habitat conservation planning at the landscape scale. This 
knowledge is germane to restoration efforts for the New England cottontail 
(Sylvilagus transitionalis), an early successional habitat specialist of conservation 
concern. We evaluated local population structure and measures of genetic 
diversity of a geographically isolated population of cottontails in the northeast- 
ern United States. We also conducted a multiscale landscape genetic analysis, in 
which we assessed genetic discontinuities relative to the landscape and devel- 
oped several resistance models to test hypotheses about landscape features that 
promote or inhibit cottontail dispersal within and across the local populations. 
Bayesian clustering identified four genetically distinct populations, with very lit- 
tle migration among them, and additional substructure within one of those 
populations. These populations had private alleles, low genetic diversity, criti- 
cally low effective population sizes (3.2-36.7), and evidence of recent genetic 
bottlenecks. Major highways and a river were found to limit cottontail dispersal 
and to separate populations. The habitat along roadsides, railroad beds, and 
utility corridors, on the other hand, was found to facilitate cottontail movement 
among patches. The relative importance of dispersal barriers and facilitators on 
gene flow varied among populations in relation to landscape composition, 
demonstrating the complexity and context dependency of factors influencing 
gene flow and highlighting the importance of replication and scale in landscape 
genetic studies. Our findings provide information for the design of restoration 
landscapes for the New England cottontail and also highlight the dual influence 
of roads, as both barriers and facilitators of dispersal for an early successional 
habitat specialist in a fragmented landscape. 



Introduction 

Understanding how landscape features influence the con- 
nectivity and genetic variation of natural populations is 
of central importance in ecology, evolution, and conserva- 
tion biology. Connectivity remains one of the most diffi- 
cult parameters to measure, yet it is a critical issue 
to address in landscape conservation (Tischendorf and 



Fahrig 2000; Lindenmayer et al. 2008). From a species' 
perspective, connectivity is a function of the ability of an 
individual to disperse through the landscape. Characteris- 
tics of habitat patches and the intervening landscape 
matrix can either facilitate or impede dispersal success 
(e.g., Perez-Espona et al. 2008). Because landscapes are 
spatially heterogeneous, and increasingly so as a result of 
human modifications, it is important to understand how 
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landscape features affect animal movement and subse- 
quent population processes. 

Landscape influences on dispersal are determined by 
species-specific characteristics, including the organism's 
vagility and specific habitat requirements for dispersal. 
These factors determine the scale and extent to which 
specific landscape features influence population connec- 
tivity. For example, broadscale dispersal barriers may 
derive from natural landforms that are impassable, such 
as mountain ranges (Zalewski et al. 2009) or ocean 
trenches (Cunningham et al. 2009), and serve to com- 
pletely separate populations. Local-scale or partial barriers 
are often formed by smaller landscape elements, such as 
roads (Coulon et al. 2006) or rivers (Frantz et al. 2010). 
The effects of these features can vary widely among spe- 
cies. Rivers may completely isolate populations of small 
mammals (Chambers and Garant 2010), but may be more 
permeable, at least under some circumstances, to larger 
mammals (Perez-Espona et al. 2008; Cullingham et al. 
2009) or even provide habitats conducive to dispersal 
(e.g., riparian corridors; Lowe and McPeek 2012). Simi- 
larly, roads, despite their recognized negative effects as 
barriers (Forman et al. 2003), may serve as movement 
corridors for some species for which associated habitat is 
conducive to dispersal (Crooks and Sanjayan 2006; Bis- 
sonette and Rosa 2009; Laurence et al. 2013). Linear land- 
scape features may have complex influences on dispersal 
even within a single species, acting as both barriers and 
facilitators of dispersal. Anthropogenic changes in land 
cover can have further consequences for connectivity, as 
habitat loss and fragmentation can impede dispersal if the 
intervening matrix is prohibitive to a species' movement 
(e.g., Keyghobadi et al. 1999; Dixon et al. 2007). These 
consequences are more pronounced for species with high 
habitat specificity (Rothermel and Semlitsch 2002). 

Disruption of habitat connectivity typically leads to 
genetic structuring among individuals, as a result of isola- 
tion (Segelbacher et al. 2003) and/or physical barriers to 
dispersal and concomitant gene flow (McRae et al. 2005). 
Reduced genetic exchange (i.e., fewer dispersing and sub- 
sequently reproducing individuals) among populations 
results in the gradual genetic divergence of populations 
through genetic drift and local adaptation (Willi et al. 
2007) or, in the extreme, leads to population extinction 
(Bond et al. 2006). The consequences of reduced connec- 
tivity are especially relevant for species of conservation 
concern, which often exist in small, isolated patches and 
have limited dispersal and small effective population sizes 
(Ewers and Didham 2006). Small populations are more 
susceptible to stochastic events, as well as a loss of genetic 
diversity, which limits the population's ability to cope 
with environmental change (Templeton et al. 1990, 2001). 
In such cases, it is important to identify gene flow 



barriers that can be mitigated to increase effective dis- 
persal. Improving connectivity helps maintain genetic 
diversity and increases effective population sizes, thereby 
strengthening the probability of population persistence 
(Newman and Pilson 1997; Frankham 2005; Bailey 2007). 
Additionally, recognizing landscape features that facilitate 
dispersal is necessary for species' recovery, so that those 
features can be maintained and replicated in habitat res- 
toration efforts to increase connectivity and augment gene 
flow where needed. 

Issues of connectivity are germane for organisms that 
rely on early successional and shrubland habitats. These 
ephemeral habitats occur in a landscape mosaic of habitats 
in varying successional stages, many of which are inhospi- 
table to early successional habitat specialists. Although pat- 
chy by nature, the spatial configuration (abundance, patch 
size, and distribution) of early successional habitats has 
been modified by a loss of natural disturbance regimes, 
land use change, and anthropogenic landscape modifica- 
tions. These habitats are on the decline in the northeastern 
United States, along with many species that rely on them 
(Brooks 2003; Litvaitis 2003; Lorimer and White 2003; 
Sauer et al. 2011). Consequently, early successional habitat 
specialists may face the consequences of habitat loss and 
fragmentation, including population isolation and decline, 
and concomitant reduction in genetic variation (Andren 
1994; Fahrig 2003; Keyghobadi 2007). 

One of the many shrubland obligate species of high 
conservation priority in the northeastern United States is 
the New England cottontail (Sylvilagus transitionalis; 
Fig. 1), which requires dense, brushy vegetation for food 
and escape cover (Litvaitis et al. 2003). Widespread habi- 
tat loss has resulted in rapid population decline for this 
species, and it now occupies less than 14% of its historical 




Figure 1. Young New England cottontail in old-field habitat. USFWS 
Photograph by Kelly Boland. 
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range (all New England states and eastern New York; Lit- 
vaitis et al. 2006). As a result, the New England cottontail 
is listed as endangered in Maine (MDIFW 2007) and New 
Hampshire (NHFG 2008), and it is a candidate for federal 
listing under the Endangered Species Act (USFWS 2006, 
2012). Remnant populations of New England cottontail 
currently occur in five geographically (Litvaitis et al. 
2006) and genetically (Fenderson et al. 2011) isolated 
regions: (1) southern Maine and Seacoast (southeastern) 
New Hampshire; (2) central New Hampshire; (3) Cape 
Cod, Massachusetts; (4) eastern Connecticut and Rhode 
Island; and (5) western Connecticut, western Massachu- 
setts, and eastern New York. The current population 
structure is a result of recent habitat fragmentation 
(within the last several decades) and genetic stochasticity, 
as the populations have experienced genetic drift in isola- 
tion (Fenderson et al. 2011). 

Given the lack of gene flow among the remaining pop- 
ulations of New England cottontails, conservation efforts 
must begin within each of these regions to ensure con- 
nectivity, stability, and population persistence on a local 
scale. New England cottontails in southern Maine and 
the Seacoast region of New Hampshire are in immediate 
need of restoration management. This region is at the 
northern extent of the species' range and is experiencing 
ongoing decline, with an estimated 50% reduction in 
effective population size occurring within the past two 
decades (Fenderson 2010) and reduced genetic diversity 
relative to other remnant populations (Fenderson et al. 
2011). A census population size of roughly 300 individu- 
als has been estimated to occur in southern Maine (Litva- 
itis and Jakubas 2004), and an effective population size of 
75-150 has been estimated for the Maine and New 
Hampshire region (Fenderson et al. 2011). Extensive hab- 
itat loss and fragmentation have reduced the availability 
of suitable (thicket) habitat in this region, such that fewer 
and smaller habitat patches exist, separated by increas- 
ingly large geographic distances. Remaining habitat 
patches are typically small (2-35 ha, with most <5 ha) 
and fragmented by development and inhospitable habitat. 
Recovery of the New England cottontail in the Maine 
and Seacoast New Hampshire region will require increas- 
ing available suitable habitat to support patch occupancy, 
as well as increasing connectivity among remaining 
patches. These efforts require an understanding of current 
landscape influences on gene flow. 

The objectives of our study were threefold: (1) to assess 
local population genetic structure and diversity of New 
England cottontails in southern Maine and coastal New 
Hampshire; (2) to identify landscape features that are 
influential in structuring populations through promoting 
or inhibiting connectivity within and among these popu- 
lations; and (3) to test hypotheses about the influence of 



landscape features (identified in #2) on gene flow. Specifi- 
cally, we evaluated the effects of geographic distance, 
roads, waterbodies, and linear features comprised of early 
successional habitat, such as utility lines and roadsides, 
on gene flow. We expected to find fine-scale population 
structure resulting from the separation of populations by 
fragmentation and/or dispersal barriers. We hypothesized 
that landscape features have a stronger influence on 
genetic variation within and among populations than geo- 
graphic distance alone. We also predicted that roads and 
waterbodies would function as dispersal barriers, while 
linear shrubby habitat features (railroads, powerline 
rights-of-way, and roadsides) would facilitate gene flow. 
Our results provide key information for the design of res- 
toration landscapes that enhance connectivity for the New 
England cottontail and thereby likely also benefit other 
species that rely on early successional habitats. Our find- 
ings also illustrate the complexity of natural and anthro- 
pogenic factors influencing gene flow of a habitat 
specialist in a fragmented landscape. 

Methods 

Study design and sample collection 

We conducted systematic fecal pellet surveys across the 
recently occupied range of New England cottontails in 
southern Maine and Seacoast New Hampshire during the 
winters of 2007/2008 and 2008/2009. Surveying in the 
winter increases detectability due to the presence of tracks 
in the snow and the increased visibility of pellets (Bru- 
baker et al. 2014). Cold temperatures promote preserva- 
tion and yield of DNA in fecal pellets (Kovach et al. 
2003). Additionally, winter sampling occurs after juveniles 
from the previous summer have dispersed (Chapman and 
Ceballos 1990) and prior to parturition of the first litter 
of the year (Chapman 1975). Sampling is thereby limited 
to postdispersal adults, and inadvertent sampling of 
highly related litter groups is avoided. The sampling 
of kin groups is further precluded by the solitary nature 
of New England cottontails (Litvaitis et al. 2008). 

Sampling scheme and scale are important consider- 
ations in planning a landscape genetics study, and they 
can influence the conclusions reached (Anderson et al. 
2010; Segelbacher et al. 2010). An ideal sampling scheme 
should incorporate the range of spatial and genetic vari- 
ability by sampling a relatively fine grain size (with 
respect to an organism's dispersal distance) across a rela- 
tively large geographic area (Storfer et al. 2007; Schwartz 
and McKelvey 2009). For our objective of identifying 
landscape influences on a fine scale, a continuously dis- 
tributed sampling scheme in areas of occupancy was 
appropriate (Storfer et al. 2007; Schwartz and McKelvey 
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2009). We surveyed 191 patches in 2007-2008 to deter- 
mine occupancy and used these occupancy results as pilot 
data to plan sampling in the subsequent field season 
(Fenderson 2010). 

Our sampling design in 2008-2009 was intended to 
obtain representative genotypes distributed continuously 
across the occupied landscape using a hierarchical system- 
atic grid pattern. To optimize search effort and the number 
of unique individuals sampled, sampling was conducted 
using finer grains (1-2 km) in areas of recent occurrence 
and coarser grains (4-8 km) as the likelihood of encoun- 
tering a New England cottontail decreased. Surveys were 
centered around grid points where we searched up to three 
suitable (densely shrubby) habitat patches within an 
approximate 1 km radius around each grid point if 



possible, although not all grid points had nearby suitable 
habitat (Fenderson 2010; Fig. 2). Within each occupied 
patch, we collected samples consisting of up to 10 pellets 
from a single pile or set of tracks, assumed to be from a 
single individual. Where possible, multiple samples were 
collected per patch, separated by at least 50 m, to maxi- 
mize the number of individuals sampled. This was the 
most exhaustive sampling effort in this area to date and 
likely documented nearly all currently occupied New Eng- 
land cottontail patches in Maine and Seacoast New Hamp- 
shire. All pellets were stored at — 20°C until analyzed. Also 
included in our dataset were three pellet samples collected 
in the winter of 2006/2007, seven opportunistically col- 
lected predator-kill or road-kill tissue samples, and blood 
samples from 19 animals trapped for relocation in 2010. 




Figure 2. Sampling scheme for surveys of 
New England cottontail fecal pellets during the 
winter of 2008/2009 and all patches searched 
during both field seasons (2007/2008 and 
2008/2009). Stars show the grid points used to 
center surveys (see text for explanation). 
Circles indicate all patches that were surveyed 
for this study. Red circles identify New England 
cottontail samples collected and yellow circles 
identify patches that were occupied by New 
England cottontails in 2000-2003, but no 
longer occupied in 2007/2008 or 2008/2009. 
Gray circles depict all of the remaining patches 
that were searched but were not occupied by 
cottontails. 
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Microsatellite genotyping 

DNA was extracted from one pellet per sample using the 
QIAamp® DNA Stool Mini Kit (Qiagen, Valencia, CA), 
with minor modifications of the manufacturer's instruc- 
tions as described in Kovach et al. (2003). As New Eng- 
land cottontails are sympatric with the eastern cottontail 
(Sylvilagus floridanus) and the snowshoe hare (Lepus 
americanus) in portions of their range, the species of ori- 
gin of the pellets was determined using a combination of 
two diagnostic RFLP tests of the mitochondrial DNA, 
using the restriction enzymes Bfal (Litvaitis and Litvaitis 
1996) and Nlalll (Kovach et al. 2003), following Fender- 
son et al. (2011) and Kilpatrick et al. (2013). 

New England cottontail samples were amplified with 
fluorescent dye-labeled primers and multiplex protocols 
at 11 microsatellite loci (Table 1) in a two-tiered 
approach. First, we used eight loci previously found to be 
polymorphic in cottontails in this study area, including a 
SRY microsatellite for sex determination (Fenderson et al. 
2011), to screen unique individuals from replicate samples 
collected within a patch. These eight-locus genotypes were 
compared in Dropout (McKelvey and Schwartz 2005) to 
identify unique individuals (P ID sws f° r the seven autoso- 
mal loci = 1.71 1e-2; Pidsibs including the SRY 
locus = 1.135E-2). Samples from unique individuals were 
then genotyped at three additional loci determined to be 
polymorphic in this study (Sfl8, Sflll, and SfllS; Berkman 
et al. 2009). PCR products were electrophoresed on an 
ABI 3130 automated DNA sequencer (Applied Biosys- 
tems, Foster City, CA). Genotypes were manually scored 
using Peak Scanner 1.0 software (Applied Biosystems), 
and alleles were binned in reference to a positive control 
with the program Allelogram 2.2 (available at http:// 
code.google.eom/p/allelogram/), to ensure consistency of 
allele calls across multiple electrophoretic runs. 

To address issues of genotyping error, PCR amplifica- 
tion and electrophoresis were replicated at least three 
times for each sample at the eight initial loci and the 



three additional loci for unique samples, until a consensus 
genotype was reached. We required alleles to amplify at 
least twice for an individual to be scored as a heterozy- 
gote at a locus. Following Frantz et al. (2003), if this rule 
was not met with the initial three replicate PCRs, we 
repeated amplifications in a stepwise fashion, for up to 
seven replicates, until each allele was observed at least 
twice. If the DNA sample was exhausted before all repli- 
cate genotypes could be obtained, we still retained a 
genotype at a locus if it successfully amplified twice and 
an identical genotype was obtained each time (only 5% of 
the 8-locus consensus genotypes and 5% of the final 11- 
locus dataset were based on two amplifications). Samples 
missing data at four or more loci were excluded from 
analyses. To ensure unique individual identity of geno- 
types, we reinspected the raw genotype peaks of all pair- 
wise samples that differed by <3 loci. Genotyping error 
was assessed by manually comparing each replicate geno- 
type to the consensus (Taberlet et al. 1996) and calculat- 
ing total error rates following Pompanon et al. (2005). 
We used Microchecker (Van Oosterhout et al. 2004) to 
test for null alleles with the Brookfield 1 estimator 
(Brookfield 1996) and Genepop (Raymond and Rousset 
1995) to test for Hardy-Weinberg equilibrium (HWE) 
and gametic disequilibrium. 

Population genetic structure, diversity, and 
effective size 

We assessed population genetic structure using two 
Bayesian clustering methods: Structure (Pritchard et al. 
2000), a program that delineates genetically similar indi- 
viduals based solely on the genetic data and not on a pri- 
ori population definitions, and Tess (Chen et al. 2007), a 
similar program that also incorporates sampling locations 
to help define genetic units. Structure was run 20 times 
at each K (the number of putative genetic populations) 
from 1-7 using a burn-in of 100,000 and run length of 
500,000, with the no admixture and independent allele 



Table 1. Multiplex PCR conditions for microsatellite loci used in this study. 



Locus/Multiplex 


Individual Primer Cone. 


MgCI 2 


dNTP Cone. 


Annealing Temp. (°C) 


Cycles 


Sat12/Lsa1 


0.20 /im 


3.00 mmol/L 


0.20 mmol/L 


55 


35 


Sat13 


0.33 /im 


1.50 mmol/L 


0.20 mmol/L 


55 


38 


Sat3 


0.33 /im 


2.50 mmol/L 


0.20 mmol/L 


58 


40 


INRA016/Y 


0.53 /im 


1.50 mmol/L 


0.25 mmol/L 


60 


35 


SfIS 1 


0.40 /im 


1.00 mmol/L 


0.20 mmol/L 


58 


30 


sfn 1/sfn 5 1 


0.24 /im 


1.00 mmol/L 


0.20 mmol/L 


58 


40 


Sol44 1 


0.20 /im 


2.50 mmol/L 


0.20 mmol/L 


62 


40 


Sol03 1 


0.20 /im 


3.00 mmol/L 


0.20 mmol/L 


52 


40 



All PCRs were in 15 /(L reactions, except as noted. All reactions used 1X BSA, 1X buffer, and 0.75 U Taq polymerase. 
1 12.5 /iL reaction. 
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frequencies model (Pritchard et al. 2000). Tess was also 
run 20 times at each K from 2-7 for 600,000 iterations, 
including a burn-in period of 100,000 sweeps. We used 
the no admixture model and a spatial interaction parame- 
ter of 0.6 (Durand et al. 2009). For each analysis, the 
optimal K was determined from the plateau of the aver- 
age lnPR(X|^C) (Structure) or the average deviance infor- 
mation criterion (DIC) of each K plotted against K 
(Tess), and from evaluation of the bar plots. Individual 
population memberships were averaged using Clumpp 
1.1.2 (Jakobsson and Rosenberg 2007) and visualized with 
Distruct 1.1 (Rosenberg 2004). Based on these results 
and geographic proximity, individuals were grouped into 
genetic clusters for the remaining analyses. 

Gene flow among populations was assessed using 
assignment tests performed in Structure, using the 
genetic clusters as prior population information and the 
same burn-in and run length as above. Migration rate 
was set at 0.05, and we tested for migrant ancestry up to 
two generations back. Assignment tests were also per- 
formed in Geneclass 2.0 (Piry et al. 2004), with the Rann- 
ala and Mountain (1997) Bayesian method and a 
threshold of 0.05. We also conducted migrant detection 
using the L_home/L_max criterion with Monte Carlo re- 
sampling (Paetkau et al. 2004) and a threshold of 0.01. 
We estimated genetic differentiation among all pairs of 
population clusters using F ST , as implemented in FSTAT 
(Goudet 1995). False discovery rate (FDR) control (Benja- 
mini and Hochberg 2000) was applied to assess signifi- 
cance for multiple comparisons using the Excel 
spreadsheet Tabulator (Verhoeven et al. 2005). 

To evaluate genetic diversity within each population 
cluster, we used Genepop to calculate the number of pri- 
vate alleles in each population. Heterozygosities, allelic 
richness, and F IS were calculated in Fstat, and private 
allelic richness was calculated in Hp-rare (Kalinowski 
2005). For multiple comparisons, we implemented FDR 
control, as above. Effective population size (N e ) of each 
genetic cluster was calculated with two single sample 
methods: a linkage disequilibrium method performed in 
LDNe (Waples 2006) and a Bayesian method performed 
in ONesAMP (Tallmon et al. 2008). We also tested for evi- 
dence of a genetic bottleneck in each of the clusters 
using several approaches. We used the M-ratio method 
with a 6 of 1, assuming a historical effective population 
size of 500, and the average parameter values identified 
by Garza and Williamson (2001) of a mean step size of 
2.8 and the percentage mutations larger than single step 
of 0.12. We also utilized the one-tailed Wilcoxon signed- 
rank test for heterozygosity excess and tested for allelic 
mode-shift (Luikart et al. 1998) with Bottleneck (Piry 
et al. 1999). We set the variance to 12 and the percent- 
age of single stepwise mutations to 0.88 (for consistency 



with the parameters used for M-ratio) and ran 1000 iter- 
ations. 

Landscape influence on gene flow 

Results of the Bayesian clustering analyses and additional 
preliminary analyses of Fenderson (2010) provided some 
insight into landscape features that may be influencing 
cottontail gene flow. Based on these results, we hypothe- 
sized that roads and waterbodies were limiting dispersal. 
Additionally, we hypothesized that several linear land- 
scape features, including powerline rights-of-way, railroad 
edges, and roadsides, were conducive to dispersal because 
they tend to be comprised of shrub habitat (Tash and Lit- 
vaitis 2007). We tested our hypotheses more explicitly by 
developing several raster cost-distance models, to calculate 
pairwise individual effective geographic distances using 
the Cost Distance tool from the Landscape Genetics Tool- 
box in ArcGIS 10 (Etherington 2011). 

We developed models for three feature classes: (1) A 
roads model tested the hypothesized barrier effects of 
roads by assigning them elevated costs relative to the 
background landscape matrix (all of the nonfeature cells 
that were assumed to have equal influence on cottontail 
dispersal). The variables of interest were the six classes of 
road, defined by traffic volume, and hypothesized to have 
increasing barrier effect with increasing traffic volume. (2) 
The surficial water model similarly evaluated the dis- 
persal-limiting effects of waterbodies (including rivers, 
streams, lakes, ponds, and coastal inlets). Variables 
included the six Strahler stream order classes, hypothe- 
sized to have increasing barrier effect with increasing 
stream width; all other waterbodies without stream order 
information were grouped with the lowest stream order 
category. (3) The facilitators model assessed the effects of 
linear strips of shrubby habitat as facilitators of dispersal, 
by assigning these features reduced costs relative to the 
background matrix. Facilitator variable costs were based 
on expert opinion (i.e., biologists familiar with the study 
organism and its habitat preferences, based on extensive 
field experience), which considered railroads, class 3-6 
road edges, class 1-2 road edges, and powerlines, in 
increasing order of hypothesized facilitating effect. In 
addition, we tested a null model of straight Euclidean dis- 
tance. GIS data sources and additional details of method- 
ology can be found in Appendix 1. 

For each of the three feature classes (roads, surficial 
water, and facilitators), we developed four types of cost 
models - binary, linear, exponential, and logarithmic - to 
evaluate the relative dispersal cost of the variables in each 
feature class. In the binary model, all feature variable 
costs were equal, but higher (for the roads and surficial 
water barrier models) or lower (for the facilitator model) 
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than the background cost. In the linear, exponential, and 
logarithmic models, the feature variable costs were ranked 
according to their hypothesized influence on dispersal 
(described above). Costs varied relative to one another in 
a linear, exponential, and logarithmic relationship, respec- 
tively, with increasingly higher costs given in the two 
types of barrier models and reduced costs in facilitator 
models on the variables hypothesized to have greater 
influence on dispersal (see Table 2 for model costs). 
Analyses were conducted for the study area as a whole as 
well as within each of the genetically distinct populations, 
as defined by the Bayesian clustering analyses. For the 
within-population analyses, we excluded the Jetport (see 

Table 2. Costs used for feature variables in raster cost-distance analy- 
sis of 12 landscape models. 

Models Binary Linear Exponential Logarithmic 



Roads 1 



Background 


1 


1 


1 


1 


Trail 


10 


2 


4 


10 


(Road Class 6) 










Unimproved 


10 


3 


9 


100 


(Road Class 5) 










Improved 


10 


4 


16 


1000 


(Road Class 4) 










Secondary 


10 


5 


25 


10000 


(Road Class 3) 










Primary 


10 


6 


36 


1 00000 


(Road Class 2) 










Interstate 


10 


7 


49 


1 000000 


(Road Class 1) 










Surficial water 2 










Background 


1 


1 


1 


1 


Stream Order 1 


10 


2 


4 


10 


Stream Order 2 


10 


3 


9 


100 


Stream Order 3 


10 


4 


16 


1000 


Stream Order 4 


10 


5 


25 


10000 


Stream Order 5 


10 


6 


36 


1 00000 


Stream Order 6 


10 


7 


49 


1 000000 


Facilitators 3 










Powerlines 


1 


1 


1 


1 


Road edges 


1 


2 


4 


10 


(Classes 1-2) 










Road edges 


1 


3 


9 


100 


(Classes 3-6) 










Railroads 


1 


4 


16 


1000 


Background 


10 


5 


25 


10000 



1 Relative costs were assigned according to the road classes. 
2 Strahler stream order class was joined to National Hydrography Data- 
set waterbody and area files based on spatial location to take into 
account drainage as well as the size of the waterbody. This was used 
to assign relative costs, and all waterbodies without stream order 
information were assigned to the lowest stream order category (e.g., 
assigned a cost of "2" in the linear model). 

3 Each facilitator variable was ranked by expert opinion according to 
its presumed utility in facilitating cottontail dispersal. 
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clustering results below), due to the small number of 
individuals sampled from a small geographic area with 
few intervening landscape features of interest. 

To evaluate the relationship between gene flow and the 
three landscape feature classes, we first used Mantel tests 
(Mantel 1967) to correlate the resulting cost distances 
with pairwise individual Rousset's a r (Rousset 2000) 
genetic distances calculated in Spagedi (Hardy and Veke- 
mans 2002). We log 10 -transformed the pairwise Euclidean 
distances and used the effective geographic distances un- 
transformed. This analysis was conducted with the ecodist 
package (Goslee and Urban 2007) in R statistical software 
(R Development Core Team 2006). We ran 10,000 per- 
mutations using the nonparametric Spearman correlation. 
Significance was assessed following Bonferroni correction 
for multiple comparisons (adjusted P-value = 0.0038 for 
a = 0.05). We compared the Mantel r values for the four 
cost models for each of the three feature classes to evalu- 
ate which cost model exhibited the best linear relationship 
with genetic distance within each population. 

Although Mantel tests are an appropriate method for 
comparing data consisting of distance measures (Legendre 
and Fortin 2010), their use in landscape genetic applica- 
tions has come under recent scrutiny (Balkenhol et al. 
2009; Guillot and Rousset 2013; Graves et al. 2013). There- 
fore, we used another complementary approach - multiple 
regression on distance matrices (MRM; Manly 1986; Lich- 
stein 2007) - to evaluate the relative importance of the 
three landscape feature classes (roads, surficial water, and 
facilitators) and Euclidean distance in influencing New 
England cottontail gene flow. MRM, like Mantel tests, can 
be used with nonindependent, pairwise genetic distance 
data. However, it provides the advantage of examining the 
influence of all input matrices simultaneously and deter- 
mining the statistical significance and relative importance 
of each variable of interest (Lichstein 2007). Balkenhol 
et al. (2009) found that, in simulations, MRM performed 
better than Mantel tests, with a good balance between type 
I error and power. MRM tests were conducted using the 
cost model (binary, linear, exponential, and logarithmic) 
with the highest Mantel r for each feature class, using 
10,000 permutations and the Spearman correlation with 
the ecodist package in R. The relative influence of each fea- 
ture class on genetic variation was further elucidated using 
the hierarchical partitioning method of Chevan and Suth- 
erland (1991). This test was conducted in the hier.part 
package for R (Walsh and MacNally 2008) using the R 2 
values generated in the MRM analyses. 

Results 

Of 610 collected samples, 376 were identified as New 
England cottontail, and these samples originated from 
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only 54 of 461 searched patches. These survey results 
revealed a significant range contraction in comparison 
with the most recent surveys of Litvaitis et al. (2006) con- 
ducted during 2000-2003 (Fig. 3). Of the 376 New Eng- 
land cottontail samples, 335 samples yielded sufficiently 
complete genotypes; 157 of those were determined to be 
unique individuals. Average raw genotyping error rates 
(our estimated genotyping error per single PCR replicate) 
across loci were 0.084 per genotype and 0.043 per allele 
(Table 3). 

Population genetic structure, diversity, and 
effective size 

The Bayesian clustering methods detected hierarchical 
population genetic structure (Figs 4 and 5). For K = 3, 
the bar graphs in both Structure and Tess identified sup- 
port for three differentiated populations: (1) a large clus- 
ter of individuals from Cape Elizabeth (Cape Elizabeth; 
Fig. 3); (2) the individuals sampled at the Portland Inter- 
national Jetport (Jetport) together with those sampled on 



Figure 3. New England cottontail locations (points) and genetic 
clusters (circles) in southern Maine and New Hampshire. The dotted 
line indicates a partial barrier, consistent with the Salmon Falls/ 
Piscataqua River, which further subdivides the Kittery West population 
into Kittery West-Maine (KW-ME) and Kittery West-New Hampshire 
(KW-NH). Inset: Estimated New England cottontail maximum range 
extents in this region circa 1960, 2003, and 2009 (this study), based 
on field surveys and historical reports. 



Table 3. Raw genotyping error rates 1 (total allele scoring mismatches 
as compared to consensus genotype) for noninvasive New England 
cottontail fecal pellet samples at 10 autosomal and one Y-chromo- 
somal microsatellite loci. 



Locus 


Per-locus error rate 


Per-allele error rate 


Sat 12 


0.188 


0.096 


Sol03 


0.161 


0.083 


Sol44 


0.057 


0.029 


Lsa1 


0.031 


0.016 


Sat 13 


0.151 


0.078 


Sat3 


0.031 


0.016 


INRA016 


0.070 


0.035 


Sfl11 


0.035 


0.018 


Sfl15 


0.046 


0.023 


Sfl8 


0.023 


0.012 


INRA326 (Y) 2 


0.061 


0.032 



1 Error rates were calculated following Pompanon et al. (2005) equa- 
tions 1 (e a = m a /2nt; per-allele error rate) and 2 (ei = m/nt; per-locus 
error rate), where m represents the number of allelic (or genotypic) 
mismatches relative to the consensus genotype, n is the number of 
individual single-locus genotypes, and t is the number of replications. 
2 Error rate for the SRY locus (INRA326) was based upon samples that 
produced an amplified product in more than one PCR run. 

the western side of 1-95, including Seacoast New Hamp- 
shire (Kittery West); and 3) all of the individuals sampled 
from east of 1-95, as well as the individuals from a patch 
that directly abutted the interstate on the western side 
(Kittery East). Structure bar graphs seemed to best sup- 
port K = 3 as above; however, the LnPD began to plateau 
at K = 5, suggesting the potential for finer-scale structure 
(Fig. 4). For Tess, at K = 4 the DIC values showed a 
slight plateau and the bar graphs stabilized, with the 
Jetport differentiated into a separate cluster from Kittery 
West, albeit with some mixed ancestry (Fig. 5). Tess 
hard-clustering results indicated a subdivision in Kittery 
West that seemed to approximate the geopolitical bound- 
ary between Maine and New Hampshire that is formed 
by the Salmon Falls and Piscataqua Rivers (data not 
shown). Further, the individuals with shared ancestry 
between the Jetport and Kittery West were sampled 
northeast of the river in Maine. Given the large geo- 
graphic distance (48 km) between the Jetport and the 
closest patch of individuals in Kittery West today, we 
considered these findings to reflect a historical connection 
between these populations and determined it more bio- 
logically meaningful to view them as separate clusters 
with respect to current occupancy patterns. That is, 
despite evidence for their recent connectivity, cottontails 
in these populations are functionally separate populations 
today with genetic drift acting independently within each 
population. Lack of current migrants detected between 
these two populations (see assignment test results below) 
further supported this separation. 
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Figure 4. Determination of the number of New England cottontail 
genetic populations (K) in Maine and Seacoast New Hampshire based 
on Bayesian clustering results from Structure and Tess. For Structure, 
the number of putative populations is frequently determined by the 
highest average LnPD or where it begins to plateau. For Tess, the 
number of putative populations is also determined where the average 
deviance information criterion (DIC) begins to plateau and/or the K at 
which the Q-plots stabilize (Figure 5). Results are shown for the entire 
study area (A) and for Kittery West (B). See Results for our 
interpretation of the number of genetic clusters. 



To test for finer-scale structure, we conducted an addi- 
tional analysis with both Structure and Tess on just the 
Kittery West individuals, using the same parameters as 
before, up to K = 4. We found support for additional 
substructure within Kittery West, comprising two popula- 
tions: the individuals sampled in Maine, northeast of the 
river (KW-ME), and those individuals sampled from Sea- 
coast New Hampshire, west of the river (KW-NH; Fig. 3). 
Combining inference from the above analyses, we con- 
cluded there are four major genetic clusters of New Eng- 
land cottontails (Cape Elizabeth, Jetport, Kittery East, and 
Kittery West), with weaker substructure within Kittery 



West comprising two subpopulations (Fig. 3). Based on 
these results, individuals were grouped for the remaining 
analyses according to the dominant genetic cluster assign- 
ment of its sampling location. Downstream analyses were 
conducted both for K = 4 (considering Kittery West as 
one population) and for K = 5 (keeping separate the two 
subdivisions in Kittery West) where we deemed it rele- 
vant. 

Microchecker analyses found three loci with null allele 
frequencies >10% in at least one population: Sol03 in Cape 
Elizabeth (11.2%), Sflll in Kittery East (14%), and Sfl8 in 
Cape Elizabeth and Kittery West (17.3% and 18.6%, 
respectively). No null alleles were detected in KW-ME or 
KW-NH. Sol03 and Sfl8 were out of Hardy-Weinberg 
equilibrium in the Cape Elizabeth population; all remain- 
ing loci and populations exhibited HWE, lending support 
to the genetic cluster designations. As null alleles have 
minor impacts on genetic distance measures (Chapuis and 
Estoup 2007), we retained these loci for downstream 
analyses. In the Cape Elizabeth population, Sol03 and Sfl8 
also exhibited gametic disequilibrium, likely due to the 
null alleles, while Satl3/INRA16 and Satl3/Sflll were not 
in equilibrium in the Kittery East population. The latter 
effect is likely due to the null alleles in Sfll 1 and a possible 
Wahlund effect in Kittery East, which includes a geograph- 
ically isolated group of cottontails in the Wells area. Link- 
age disequilibrium is often found in small populations, 
especially as a result of recent isolation and subdivision 
(Frankham et al. 2002; Zartman et al. 2006). 

Genetic diversity measures were similar for each of the 
genetic clusters (Table 4). Private alleles were identified in 
each cluster except KW-NH. F IS values were significantly 
higher than zero for Sol03, Sflll, and Sfl8, due to the null 
alleles, leading to significant F IS in the Cape Elizabeth, 
Kittery East, and Kittery West populations overall. When 
calculated without the three null allele loci, F IS was not 
significant for any population. All pairwise F ST values 
were significant (overall F S t = 0.127, P < 0.001), and the 
largest differences occurred in comparisons of KW-NH 
and the other populations (Table 5). 

The two Kittery West subpopulations were separated 
for assignment tests and migrant detection. The Geneclass 
assignment test assigned 87.3% of the individuals back to 
their sampled location (quality index = 82.75%), and only 
six individuals were cross-assigned to other populations 
with relatively high probability (>75%). Only two individ- 
uals were identified as migrants by both Geneclass and 
Structure. One was an individual sampled in Kittery East 
that assigned to KW-NH, and the other was sampled in 
KW-ME and assigned to KW-NH. Five other individuals 
were assumed to have admixed ancestry based on meeting 
at least two of the following criteria: (1) >50% Geneclass 
assignment probability to a cluster other than that of 
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Figure 5. Individual assignment probabilities of New England cottontails to genetic clusters determined by a) Structure and b) Tess for (A) K = 3 
and (B) K = 4. Geographic sampling locations are indicated below each pair of graphs in bold. Genetic cluster assignments are indicated above 
each graph. Subdivisions of the Kittery West population are shown in C). 



geographic origin; (2) identified as a putative first-genera- 
tion migrant with Geneclass migrant detection; (3) <50% 
Structure resident probability; and/or (4) >10% Struc- 
ture immigrant probability. Three additional individuals 
from Cape Elizabeth were cross-assigned to the Jetport 
with >85% probability in Geneclass; however, they had 
>85% resident probability in Structure. Due to the prox- 
imity of the two populations, we also considered them as 
potentially admixed. 

Effective population sizes for each cluster ranged from 
only 3.2 in the Jetport to 36.7 in Cape Elizabeth (sample 
sizes were too small to test the Kittery West subpopula- 
tions separately; Table 6). Estimates obtained by the two 
methods were significantly different for Kittery East 
(based on nonoverlapping 95% confidence intervals). The 



Cape Elizabeth population showed signs of having experi- 
enced a recent genetic bottleneck (Table 6). It exhibited 
significant heterozygosity excess by the Bottleneck 
method under the I.A.M. and T.P.M. mutation models, 
with the Wilcoxon one-tailed probability test. Kittery East 
and KW-ME also had significant heterozygosity excess 
under the I.A.M. model, and KW-ME had a shifted allelic 
mode distribution as well. The M-ratio method also 
detected a significant genetic bottleneck in KW-ME. 

Landscape influence on gene flow 

For the Mantel tests comparing our hypothesized cost 
distances with cottontail genetic distance across the 
entire study area, all of the dispersal models, except the 
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Table 4. Genetic diversity of New England cottontail loci and genetic clusters in southern Maine and New Hampshire. Alleles, allelic richness, 
observed heterozygosity (Ho), unbiased expected heterozygosity (UHe), and F B are across all samples (per-locus data) or averaged across loci (per- 
population data). Private alleles are the total number (and private allelic richness is the sample size adjusted proportion) of private alleles for all loci 
in each population. 



Locus 




Alleles 




Ho 




UHe 






Sat12 




6 




0.689 




0.758 




0.091 


Sol03 




6 




0.571 




0.750 




0.238* 


Sol44 




4 




0.588 




0.582 




-0.011 


Sat13 




6 




0.575 




0.598 




0.098 


Lid I 




3 




U.ZU I 




U.ZZ3 




U.Ujo 


Sat3 




5 




0 273 




0 339 




0 1 96 


MM r\MU I O 




2 




0 276 




0 248 




— 0 113 


Sfl1 1 




A 
H 




0 356 




0 462 




0 229* 


Jill J 




Z 




UM I u 




n asm 




U. I j4 


Sfl8 




2 




0.159 




0.385 




0.589* 






Allelic richness 








Private alleles 






Population (A/) 


Alleles 


4 pops. 


5 pops. 


Ho 


UHe 


4 pops. 


5 pops. 




Cape Elizabeth (84) 


3.2 


2.99 


2.70 


0.443 


0.484 


3 (0.25) 


3 (0.24) 


0.086*/-0.034 


Jetport (19) 


2.5 


2.49 


2.21 


0.337 


0.323 


3 (0.29) 


3 (0.15) 


-0.044/-0.045 


Kittery East (28) 


3.1 


2.99 


2.68 


0.384 


0.433 


1 (0.13) 


1 (0.18) 


0.1 1 5*/0.076 


Kittery West (26) 


3.0 


2.85 




0.379 


0.424 


2 (0.22) 




0.110*/0.037 


KW-ME (10) 


2.6 




2.52 


0.441 


0.462 




1 (0.19) 


0.047/-0.032 


KW-NH (16) 


2.8 




2.35 


0.340 


0.334 




0(0) 


-0.021/-0.021 



♦Indicates P < 0.05 after Bonferroni correction. 

'Population level F, s values are given for both the full 10-locus dataset and without the three loci with null alleles (Sol03, Sf 1 11, and Sf 18), before 
and after the forward slash, respectively. 



Table 5. Pairwise F S t among New England cottontail genetic clusters 
in Maine and Seacoast New Hampshire. 





Cape Elizabeth 


Jetport 


Kittery East 


KW-ME 


Cape Elizabeth 










Jetport 


0.087 








Kittery East 


0.112 


0.102 






KW-ME 


0.087 


0.071 


0.127 




KW-NH 


0.165 


0.231 


0.244 


0.171 



All F57- values were significant at the 5% level. 



logarithmic road model, were statistically significant. The 
significant Mantel correlations ranged from 0.1915 to 
0.2159 and were slightly stronger than the correlation 
with Euclidian distance (r M = 0.1913); however, all confi- 
dence intervals overlapped (Fig. 6). 

Within the genetic clusters, the results of the Mantel 
tests were more varied and no single cost model per- 
formed the best for all feature classes (Fig. 6). For the 
roads feature class, the linear cost models performed the 
best in each population, while the logarithmic road mod- 
els were always the least correlated with genetic distance. 
For the surficial water and facilitator models, each cost 
model performed best in at least one of the populations. 
In Cape Elizabeth, Kittery East, and Kittery West, all 



models were significant. In Cape Elizabeth, all facilitator 
models had stronger correlations with gene flow than did 
Euclidean distance, and in Kittery East, only the linear 
surficial water and linear facilitator models were slightly 
more correlated with cottontail genetic distance than was 
Euclidean distance. For Kittery West, except for the loga- 
rithmic models, all of the facilitator and surficial water 
models explained more variation in genetic distance than 
Euclidean distance alone. For the Kittery West subpopula- 
tions, very few models were significant, likely due to low 
sample size. Only the linear road model performed better 
than the Euclidean model in KW-ME, whereas in KW- 
NH, most of the facilitator models, as well as the expo- 
nential and logarithmic surficial water models, performed 
better than the Euclidean model, but only the linear facil- 
itator model was significant. 

The MRM analysis allowed for quantitation of relative 
importance of the features in each population. The full 
models (which included the cost model with the highest 
Mantel correlation for each of the three feature classes 
and the Euclidean distance model) were significant across 
all populations and within each population, except KW- 
NH, although they explained a small amount of the total 
genetic variation (0.3-10%). For the study area as a 
whole, only the road variable had a significant positive 
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Table 6. Estimated effective population sizes (mean and 95% CI of N e ) using the LDNe and ONeSAMP estimators and results of genetic bottle- 
neck tests using Bottleneck and M-ratio methods for New England cottontail genetic clusters in southern Maine and Seacoast NH. 

Bottleneck 2 

Effective population size estimates 1 Wilcoxon Test Probability 

Mean W e (95% CI) Mean N e (95% CI) Bottleneck Mc M-ratio 

Population LDNe ONeSAMP 1AM TPM SMM Mode-Shift M-ratio 3 (W e = 500) Probability 



Cape 


36.7 (22.0-67.9) 


35.0 (24.0-71.9) 


0.001* 


0.012* 


0.097 


Normal L-shaped 


0.865 


0.815 


0.21 


Elizabeth 




















Jetport 


3.2 (1.6-13.9) 


14.1 (11.1-18.6) 


0.326 


0.787 


0.898 


Normal L-shaped 


0.829 


0.798 


0.13 


Kittery East 


5.3 (2.5-10.9) 


28.3 (21.0-59.0) 


0.016* 


0.313 


0.348 


Normal L-shaped 


0.828 


0.802 


0.11 


Kittery West 4 


16.8 (6.1-123.7) 


23.9 (17.2-50.0) 


0.138 


0.615 


0.754 


Normal L-shaped 


0.814 


0.802 


0.08 


KW-ME 






0.042* 


0.313 


0.423 


Shifted Mode 


0.781 


0.785 


0.04* 


KW-NH 






0.754 


0.947 


0.958 


Normal L-shaped 


0.830 


0.797 


0.15 



*P< 0.05. 

'Effective population sizes could not be estimated for the Kittery West subpopulations due to low sample size. 

2 For the tests performed in Bottleneck, the Wilcoxon one-tailed probability of heterozygosity excess for three mutation models (IAM = infinite allele 
model; TPM = two-phase model; and SMM = stepwise mutation model) is given, as well as results of the allelic mode-shift test. 
3 The M-ratio for each genetic cluster is specified; critical M values (Mc) were calculated using N e = 500; the M-ratio probability is the probability 
that the M-ratio is significantly lower than the Mc value. 

4 Bottleneck probabilities for Kittery West are from a separate simulation as the other five populations. 



association with genetic distance, and Euclidean distance 
had a significant, but negative, association (Table 7). For 
the analyses within the genetic clusters, however, the road 
variables were not significant in any population. Facilita- 
tor features were highly significant for Cape Elizabeth and 
Kittery West, as well as for KW-NH and marginally so 
for KW-ME. Euclidean distance had a significant negative 
correlation in Kittery West and a significant positive cor- 
relation in Kittery East. Surficial water was also positively 
associated with genetic distance for Kittery West. Hierar- 
chical partitioning of the independent effects of each of 
the features showed that, across the study area, nearly half 
of the explained variance in cottontail genetic distance 
was due to the influence of roads, and the independent 
effects of geographic distance, surficial water, and facilitat- 
ing habitat were about equal (Fig. 7). Within the genetic 
clusters, facilitating features explained the greatest per- 
centage of the genetic variation for all populations except 
Kittery East, which showed a strong influence of Euclid- 
ean distance alone. 

Discussion 

Habitat loss and fragmentation can alter the genetic struc- 
ture and diversity of natural populations through a dis- 
ruption of gene flow and metapopulation processes 
(Gonzalez et al. 1998). Effects are most pronounced in 
species with strong habitat associations, for which frag- 
mentation impedes dispersal (e.g., Rothermel and Sem- 
litsch 2002). New England cottontail populations have 
been declining for decades as a result of ongoing loss and 
fragmentation of early successional habitats (Litvaitis 



et al. 2006, 2008). The results of this study suggest that 
reduced occupancy is associated with low genetic connec- 
tivity among fragmented populations of New England 
cottontails in Maine and Seacoast New Hampshire. The 
current distribution of New England cottontails in this 
region represents a substantial range contraction since 
previous surveys in 2000-2003, when cottontails were 
found as far north as Cumberland, Maine (20 km north 
of the current northernmost location in Cape Elizabeth, 
Maine), and also occupied patches in the intervening 
region between the three disjunct, currently occupied 
areas (Litvaitis et al. 2006; Figs 2 and 3). This range con- 
traction, combined with our findings of genetically iso- 
lated populations with low genetic diversity, emphasizes 
the immediacy of restoration needs for New England cot- 
tontails in Maine and New Hampshire. 

Genetic structure, diversity, and bottlenecks 

New England cottontails in Maine and Seacoast New 
Hampshire were structured into five genetically distinct 
and geographically separated populations, the boundaries 
of which coincided with major highways, urban develop- 
ment, and rivers. Although Bayesian clustering results 
indicate recent, historical connections, gene flow is cur- 
rently absent or very minimal among these populations, 
as evidenced by assignment tests and the relatively strong 
differentiation (significant F ST values) among all pairs of 
populations. The presence of private alleles in each popu- 
lation further suggests that rapid genetic drift is occurring 
in the absence of dispersal. The distances separating 
populations greatly exceed the estimated mean dispersal 



1864 



© 2014 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



L. E. Fenderson ef al. 



Gene Flow in the New England Cottontail 



Figure 6. Mantel r correlations and 95% 
confidence intervals of genetic distance and 
effective geographic distance for cost models 
testing the influence of three types of 
landscape features on gene flow in New 
England cottontails. Roads and surficial water 
models tested hypothesized barrier effects of 
these features on dispersal, and the facilitator 
models tested the hypothesized influence of 
linear conduits in promoting dispersal. For each 
feature class, four cost models were tested 
(from left to right): binary, linear, exponential, 
and logarithmic, evaluating the relative cost of 
each feature variable on dispersal. For 
comparison, correlation with Euclidean 
distance alone (a model of isolation by 
distance) is shown. Statistical significance 
(indicated by asterisks) was assessed with 
10,000 permutations and two-tailed P < 0.05 
following Bonferroni correction (corrected 
P < 0.0038). Arrows indicate the model within 
each feature class with the highest Mantel r. 
Top panel - across the Maine-New Hampshire 
study area as a whole; bottom panel - within 
each genetically distinct population. 
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distances of New England cottontails (500 m— 3 km; Lit- 
vaitis and Villafuerte 1996; Fenderson 2010) and even the 
maximal dispersal distances of other lagomorphs (12- 
17 km; Gillis and Krebs 1999; Estes-Zumpf and Rachlow 
2009; Bray et al. 2007). While long-distance dispersal 
events are important in population dynamics of small 
mammals (Diffendorfer and Slade 2002), it is unlikely 
that cottontails can disperse the current interpopulation 
distances necessary to maintain gene flow among remnant 
populations. 

Fragmentation and subsequent population isolation 
have had negative genetic consequences for New England 
cottontails in Maine and Seacoast New Hampshire. Fen- 
derson et al. (2011) found that genetic diversity, as mea- 
sured by allelic richness (Ar = 2.6-2.9) and heterozygosity 
(H 0 = 0.223-0.287), was reduced in Maine and New 
Hampshire relative to geographic areas in the core of the 
species' range (Ar = 3.4-4.0 and H 0 = 0.371-0.492 for 
eastern Connecticut/Rhode Island and western Connecti- 
cut/New York), while common ancestry estimates (F- val- 
ues) were increased (0.19 in Maine/New Hampshire 
compared to 0.09-0.12). Genetic diversity was similarly 
low across the four populations identified in this study, 
with slightly reduced allelic richness and heterozygosity in 
the Jetport relative to the other populations. The reduced 
genetic diversity in the Jetport is consistent with the ori- 
gin of these 19 individuals from a single habitat patch in 



an isolated area. Genetic diversity has been found to have 
important effects in determining population dynamics 
(Reed et al. 2007), warranting further investigations into 
the implications of low genetic diversity on individual fit- 
ness and potential inbreeding effects on cottontail popula- 
tions. Concerns about inbreeding may be germane in 
light of our subsequent research findings of high genetic 
similarity among cottontails on some small patches in 
Maine and New Hampshire (Brubaker 2012; A. Kovach 
unpubl. data). 

Further genetic consequences are evidenced by genetic 
bottleneck tests, which indicated that several of the popu- 
lations have recently experienced a demographic bottle- 
neck or possibly are currently undergoing one. Of the 
two methods we used, the Bottleneck approach is more 
sensitive for detecting recent bottlenecks (within a few 
dozen generations), while the M-ratio test is best for 
detecting more historic or longer-duration bottlenecks 
(Williamson-Nateson 2005). Our results are most consis- 
tent with recent bottleneck effects, with the exception of 
KW-ME, which showed significance with both Bottleneck 
methods and the M-ratio test. This might indicate that 
the bottleneck effects are most severe in this population, 
which is bounded on the west by the river and on the 
east by the interstate and is now effectively isolated from 
the northern population. The evidence for a recent bottle- 
neck was also strong in Cape Elizabeth, the northernmost 
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Table 7. Multiple regression on matrices (MRM) analysis of the influence of three types of landscape features - roads (Rd), surficial water (Rvr), 
and facilitating habitat (Facil) - and Euclidean distance (Euclid) on New England cottontail gene flow across and within five populations in south- 
ern Maine and New Hampshire. The full models were constructed using the cost model - binary (Bnry), linear (Lnr), exponential (Exp), and loga- 
rithmic (Log) - with the highest Mantel correlation for each feature class. P? values are given for each full model and ft values for each variable in 
the model. Significant P values (P < 0.05) are indicated in bold. 
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Figure 7. Hierarchical partitioning of the 
independent effects of Euclidean distance and 
three landscape feature class types on New 
England cottontail genetic distance across the 
study area and within each population in 
southern Maine and Seacoast New Hampshire. 
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New England cottontail population rangewide, as even 
the highly conservative stepwise mutation model 
approached significance. The lack of significance with the 
M-ratio method in the Cape Elizabeth population sup- 
ports a recent demographic decline, consistent with our 
documented range contraction in the last decade, its cur- 
rent separation from the nearest population to the south 
by 29 km, and its isolation from the north and west by 
interstates. Populations at the periphery of a species' 
geographic range often have reduced gene flow, genetic 
variation, and effective population sizes (Schwartz et al. 
2003), which are often manifest in genetic signatures of 
bottlenecks. 

Along with low genetic diversity and bottleneck signa- 
tures, we found evidence of low effective population sizes 
for cottontails in this region. Previously, Fenderson et al. 
(2011) estimated the effective population size for the 
entire Maine/New Hampshire population (including a 
small cluster of individuals in central New Hampshire) to 
be 75-150 individuals. With further analysis, we have 
found that cottontails in this region actually occur in sev- 
eral small populations with critically low effective popula- 
tion sizes of <40 individuals. The lower the effective 
population size, the greater the likelihood of negative 
genetic consequences, such as inbreeding and extinction 
through stochastic effects (Allendorf and Ryman 2002). 
Conventionally, an effective population size of at least 50 
is suggested for short-term persistence, while an effective 
size of 500 is considered necessary to maintain long-term 
evolutionary potential (Franklin 1980; Franklin and 
Frankham 1998; Jamieson and Allendorf 2012). Lag- 
omorphs, however, may require an effective size of >300 
to even persist for 40 generations (based on a census size 
of 3000 - Newmark 1987; Reed et al. 2003; assuming a 
conservative N e /N ratio of 0.1, Frankham 1995). Within 
this context, New England cottontails in the Maine and 
Seacoast New Hampshire region currently do not exist in 
populations large enough to persist into the near future. 

Lagomorphs as a taxon may be particularly vulnerable 
to extinction, likely due to their short generation times 
and large fluctuations in population size, despite their 
high growth rates (Newmark 1995). The survival advanta- 
ges for species with high growth rates persist only at large 
population sizes, and high growth rate species have higher 
extinction risk than lower growth rate species at small 
population sizes (Pimm et al. 1988; Newmark 1995). Low 
population sizes may act synergistically with poor habitat 
quality, such as that resulting from anthropogenic influ- 
ence, to further increase extinction vulnerability (Reed 
et al. 2003). Reduced genetic variation and effective pop- 
ulation size may negatively impact survival, fecundity, 
and population growth rates (Reed et al. 2007). These cir- 
cumstances are important to consider in predicting the 



future persistence of New England cottontails in Maine 
and Seacoast New Hampshire, where the limited, 
degraded, and fragmented suitable habitat, combined with 
reduced genetic diversity, likely exacerbates the vulnerabil- 
ity of these small populations. 

Anthropogenic and natural influences on 
gene flow 

Extensive movements in a fragmented landscape likely 
come at significant costs in the form of increased energy 
expenditure and high mortality risks. Even distances 
>5 km may be difficult for cottontails to overcome, as we 
previously found genetic discontinuities associated with 
this level of patch isolation (Fenderson 2010). Successful 
dispersal among disjunct patches is likely strongly depen- 
dent on the intervening landscape matrix. In this study, 
we found that three matrix features - roads, waterbodies, 
and linear conduits of thicket habitat - influenced gene 
flow of cottontails. The relative importance of each fea- 
ture type, however, was a function of the landscape 
matrix at the scale of analysis and varied by population, 
illustrating the effects of scale and landscape gradients 
(Schwartz and McKelvey 2009; Cushman and Landguth 
2010; Jaquiery et al. 2011). 

At all spatial scales, roadsides and other facilitating 
habitat features had positive effects, while roads, water- 
bodies, and geographic distance had negative effects on 
cottontail gene flow. Across the study area as a whole, 
major highways, the river, and geographic isolation subdi- 
vided cottontail populations, while within populations, 
features that facilitate dispersal between suitable habitat 
patches were important in maintaining gene flow on a 
local scale. In Cape Elizabeth, where occupied patches 
were large and proximate and the landscape matrix con- 
tained few dispersal-limiting features (few waterbodies 
and only low-volume roads), only facilitating habitat was 
important in explaining gene flow. In Kittery East, which 
is the most fragmented population, comprised of two dis- 
junct clusters of individuals separated by 20 km and 
where remnant patches are small (average patch size is 
2.1 ha compared to 3.9 ha in Kittery West and 5.4 ha in 
Cape Elizabeth; Fenderson 2010), geographic distance was 
the most important factor. In Kittery West, the dispersal- 
limiting effect of a large river dominated, with facilitating 
habitat also influential in explaining gene flow. For sub- 
populations within Kittery West, results were variable 
among the different analytical methods and no clear land- 
scape pattern emerged. This variability may have been a 
result of small sample sizes. Alternatively, at these small 
spatial scales, dispersal patterns may be more influenced 
by microsite characteristics or behavioral interactions 
between individuals. 
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Our clustering analyses suggested a barrier effect of the 
two interstates in our study area, 1-95 and 1-295, and in 
our MRM resistance modeling, the roads-as-barriers 
model explained the largest amount of genetic variation 
in the analysis across all populations. The size of the road 
is an important factor influencing dispersal, and it is 
likely that only major roads with high traffic loads are 
substantial barriers to dispersal (e.g., Frantz et al. 2010). 
Nevertheless, the logarithmic road model, which placed 
an extremely high dispersal cost on the interstate, was the 
only road model that was not significant in the Mantel 
tests across all populations, suggesting that even major 
highways are not absolute barriers. Further, the linear cost 
model, for which costs increase incrementally with road 
class, was the top road model for analyses across and 
within all populations. Accordingly, in support of an 
incomplete barrier effect of roads, we found genetic simi- 
larity of cottontails that occupied patches on either side 
of 1-95 in Kittery, adjacent to the highway. Dispersal 
between these patches may have occurred through a cul- 
vert that passes underneath the highway in this location. 
Alternately, this connectivity might be a result of an occa- 
sional individual successfully crossing the highway, as has 
been observed by radiotelemetry (J. Litvaitis, pers. obs.; 
H. Holman, New Hampshire Fish and Game, pers. 
comm.). Underpasses with shrubby riparian habitat may 
also facilitate cottontail dispersal across interstates. Such 
an underpass occurs in the vicinity of the Portland Jet- 
port and, in combination with historical occupancy of 
previously suitable habitat patches (discussed below), may 
explain the genetic connectivity observed between the Jet- 
port individuals (east of 1-95) and those in Kittery West 
(west of 1-95), despite the barrier posed by the interstate. 
Large highways have been found to restrict movement in 
other small- to mid-sized mammals, such as badgers (Me- 
les meles; Mata et al. 2008) and pygmy rabbits (Bmchyla- 
gus idahoensis; Estes-Zumpf et al. 2010), although they 
may use permeable features, including culverts and 
underpasses (Mata et al. 2008). 

Despite the widely recognized negative ecological effects 
of roads (Balkenhol and Waits 2009; Fahrig and Rytwin- 
ski 2009), they appear to have a complex effect on natural 
populations that may vary with focal species, population 
size, and road type (Clevenger et al. 2001; Forman et al. 
2003; Gauffre et al. 2008). Roads, which are often associ- 
ated with adjacent strips of herbaceous and shrubby vege- 
tation, can create and enhance habitat for some species 
(Bissonette and Rosa 2009; Fahrig and Rytwinski 2009) 
and thereby serve as movement and dispersal corridors 
rather than barriers (Crooks and Sanjayan 2006). Indeed, 
roads may enhance gene flow for some generalist and 
invasive species (Crispo et al. 2011; Laurence et al. 2013). 
Species that specialize on early successional habitats, 



including the New England cottontail, however, may be 
faced with conflicting positive and negative effects of 
roads, which may facilitate dispersal through suitable 
roadside habitat, while simultaneously increasing mortal- 
ity risk through road crossings (Tash and Litvaitis 2007). 
Given these dual facilitating and barrier effects, interstate 
highways may have an effect similar to that of drift fences 
for cottontails, which may avoid the high volume roads 
and be more likely instead to travel along them, utilizing 
the adjacent shrubby habitat to avoid crossing them (e.g., 
Forman et al. 2003; Holderegger and Di Giulio 2010). 

We were able to evaluate the roads-as-dispersal-facilita- 
tors hypothesis, in part, through our landscape resistance 
modeling, in which the facilitator model accounted for the 
facilitating effects of roadsides and other linear thicket 
conduits. The facilitator models were significant in most 
of the within-population analyses and explained a larger 
portion of the genetic variation than the barrier effects of 
roads, surficial water, or geographic distance within popu- 
lations, except for Kittery East, which showed a pro- 
nounced pattern of isolation by distance. The potential 
facilitating effects of linear shrub-lined conduits were fur- 
ther illustrated by the Bayesian clustering results of genetic 
similarity of the Jetport and Kittery West populations, 
separated by a distance of 48 km. This genetic similarity 
likely reflects recent historical connectivity. A major pow- 
erline runs parallel to and on the west side of the interstate 
between these two populations, and cottontails occupied 
habitat patches within this intervening area as recently as 
2000-2003 (Litvaitis et al. 2006; Fig. 2). The shrubby habi- 
tat along this powerline and along the interstate itself may 
have served as a north-south dispersal corridor, connect- 
ing these patches in the recent past. Although the habitat 
patches between Kittery West and the Jetport were either 
no longer suitable or unoccupied during our surveys, we 
found cottontails within and adjacent to the powerlines 
within Kittery West. While our results highlight the 
importance of linear conduits as dispersal facilitators, our 
approach did not allow us to fully evaluate the relative 
importance of the various facilitating features (roadsides, 
powerlines, and railroads), as we did not test different per- 
mutations of the relative costs for each feature. Our find- 
ings suggest that the relative importance of the facilitating 
features may depend on the landscape matrix composi- 
tion, as the best-supported facilitator cost model (highest 
Mantel r) varied by population. Teasing apart these influ- 
ences is a potentially important avenue for future research. 
Additionally, although our analyses focused on shrubland 
habitats, future studies should investigate the potential 
facilitating effects of other types of early successional habi- 
tats, such as tall grasslands and hayfields, and fully evalu- 
ate the role of other habitats and land-cover features in 
influencing cottontail dispersal. 
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Our landscape models and Bayesian clustering results 
also revealed the influence of the Salmon Falls/Piscataqua 
River as a partial barrier in the Kittery West population, 
separating individuals to its east and west. The Piscataqua 
River ranges 250-500 m wide where it empties into the 
Atlantic Ocean and is approximately 50 m wide to the 
northwest, where it becomes the Salmon Falls River. One 
of the two migrants identified by assignment tests was a 
putative disperser across the river (sampled on the east 
side in KW-ME and assigned to KW-NH on the west 
side), suggesting that cottontails do disperse, at least occa- 
sionally, across the river. Although rivers pose barriers to 
the dispersal of many small- and medium-sized mammals 
(e.g., Cullingham et al. 2009; Chambers and Garant 2010; 
Frantz et al. 2010), several species of rabbits, including 
other Sylvilagus, have been reported to swim (Chapman 
and Feldhamer 1981; Chapman and Willner 1981; Estes- 
Zumpf and Rachlow 2009). The permeability of different 
waterways may vary, however, according to their width, 
flow, winter ice cover, or surrounding landscape shape 
(Cullingham et al. 2009; Frantz et al. 2010). Small islands 
that occur in the narrower portion of this river may fur- 
ther facilitate occasional crossing by New England cotton- 
tails. These findings also bear relevance to understanding 
the distribution of nonnative eastern cottontails (Sylvila- 
gus floridanus) in this region, which extended their range 
into southern New Hampshire in the late 1960s (Jackson 
1973). Our results support the contention that the Piscat- 
aqua River has been a dispersal barrier preventing the 
spread of eastern cottontails from southern New Hamp- 
shire into Maine. Yet, these findings also raise concerns 
about the potential for eastern cottontails to cross the 
narrower portions of this river if they continue their 
spread farther northeast into New Hampshire. 

Conclusion 

This study found negative genetic consequences of frag- 
mentation and influences of landscape structure on gene 
flow for a habitat specialist. Our findings of isolated pop- 
ulations with low effective population sizes and low 
genetic diversity suggest that the New England cottontails 
in Maine and Seacoast New Hampshire are vulnerable to 
extirpation without immediate human intervention. 
Extensive habitat loss and fragmentation have reduced the 
availability of suitable thicket habitat in this region, such 
that fewer and smaller habitat patches exist, separated by 
increasingly large geographic distances. As a result, occu- 
pancy has declined, and remaining cottontails are effec- 
tively isolated into small populations, within which 
genetic drift occurs and genetic diversity is being lost in 
the absence of gene flow. Genetic data revealed historical 
connections among remnant populations, a finding that 



points toward the importance of restoring suitable habitat 
to reconnect these populations. Landscape resistance 
models also showed the importance of linear conduits of 
thicket habitats (powerlines, roadsides, railways) in sus- 
taining gene flow and the role of major highways and 
waterways in impeding dispersal. We also found evidence 
that anthropogenic connections, such as underpasses and 
possibly culverts, may be effective in facilitating dispersal 
across interstate highways. 

Management to create additional suitable habitat is 
critical for restoration of cottontail populations in this 
region. This habitat creation has been the dominant focus 
of a recent conservation initiative. The current goals of 
the conservation strategy for the New England Cottontail 
(Fuller and Tur 2012) outline targets for the size, number, 
and proximity of restored habitat patches per each desig- 
nated focal management area. Given the critically low 
effective population sizes, however, habitat creation alone 
may be an insufficient management solution and translo- 
cations may be necessary to augment existing populations. 
In addition, the creation of dispersal corridors, such as 
expanding roadside shrubby edge and potentially mitigat- 
ing highway crossings via underpasses or culvert modifi- 
cations (e.g., Dodd et al. 2004), may also be effective in 
restoring connectivity in this highly fragmented landscape. 
Our findings in this study highlight the need for consid- 
ering not only the number of hectares restored, but also 
the placement and configuration of habitat patches to 
afford gene flow within restoration landscapes. Our 
results provide a starting point for addressing the broader 
goal of designing conservation landscapes that support 
viable, functionally connected metapopulations with the 
potential to persist in the long term. Doing so will require 
establishing functional connections both within and 
among focal management areas. 

Acknowledgments 

Funding for this research was provided by the U.S. Fish 
and Wildlife Service; Maine Department of Inland Fisher- 
ies and Wildlife - through its Pitman Robertson Funds, 
State Wildlife Grant, and Department of Transportation 
Grant; the Maine Outdoor Heritage Fund; the National 
Science Foundation (Grant #BCS-126301); and the New 
Hampshire Agricultural Experiment Station. Numerous 
biologists from Rachel Carson National Wildlife Refuge 
were instrumental in assisting with sample collection, 
including David Tibbetts, Kristina Vagos, and Margaret 
Arbuthnot. Laboratory space was generously provided by 
David Berlinsky. Elisha Allen, Melanie Schroer, Grace 
Smarsh, Samantha Petren, Cynthia Sirois, and Allison 
Citro assisted with DNA extractions and species identifi- 
cations. Thanks to Jennifer Walsh for screening the Sfl 



2014 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



1869 



Gene Flow in the New England Cottontail 



L. E. Fenderson ef al. 



primers and to Niko Balkenhol and Ken Anderson for 
their helpful assistance with statistical analysis questions. 
Stephanie Coster, Jennifer Walsh, Katrina Papanastassiou, 
and two anonymous reviewers provided helpful sugges- 
tions that improved this manuscript. The findings and 
the conclusions of this article are those of the authors 
and do not necessarily reflect the view of the U.S. Fish 
and Wildlife Service. This is Scientific Contribution Num- 
ber 2455 of the New Hampshire Agricultural Experiment 
Station. 

Data Accessibility 

Microsatellite genotypes and sampling locations are 
available in Dryad (doi:10.5061/dryad.ls834) following a 
2-year embargo period. 

Genotyping protocols and cost values used in landscape 
resistance modeling are provided in manuscript tables; 
additional information on the derivation of GIS layers for 
landscape modeling is provided in the Appendix 1. 

Conflict of Interest 

None declared. 
References 

Allendorf, F. W., and N. Ryman. 2002. The role of genetics in 
population viability analysis. Pp. 50-85 in S. R. Beissinger 
and D. R. McCullough, eds. Population viability analysis. 
Univ. of Chicago Press, Chicago. 

Anderson, C. D., B. K. Epperson, M. J. E. Fortin, R. 

Holderegger, P. M. A. James, M. S. Rosenberg, et al. 2010. 
Considering spatial and temporal scale in landscape-genetic 
studies of gene flow. Mol. Ecol. 19:3565-3575. 

Andren, H. 1994. Effects of habitat fragmentation on birds and 
mammals in landscapes with different proportions of 
suitable habitat: a review. Oikos 71:355-366. 

Bailey, S. 2007. Increasing connectivity in fragmented 

landscapes: an investigation of evidence for biodiversity gain 
in woodlands. For. Ecol. Manage. 238:7-23. 

Balkenhol, N., and L. P. Waits. 2009. Molecular road ecology: 
exploring the potential of genetics for investigating 
transportation impacts on wildlife. Mol. Ecol. 18:4151-4164. 

Balkenhol, N., L. P. Waits, and R. J. Dezzani. 2009. Statistical 
approaches in landscape genetics: an evaluation of methods 
for linking landscape and genetic data. Ecography 32:818- 
830. 

Benjamini, Y., and Y. Hochberg. 2000. On the adaptive control 
of the false discovery rate in multiple testing with 
independent statistics. J. Educ. Behav. Stat. 25:60-83. 

Berkman, L. K., M. J. Saltzgiver, E. J. Heist, C. K. Nielsen, C. 
L. Roy, and P. D. Scharine. 2009. Hybridization and 
polymorphic microsatellite markers for two lagomorph 



species (Genus Syivilagus): implications for conservation. 

Conserv. Genet. Resour. 1:419-424. 
Bissonette, J. A., and S. A. Rosa. 2009. Road zone effects in 

small-mammal communities. Ecol. Soc. 14:27. 
Bond, J. E., D. A. Beamer, T. Lamb, and M. Hedin. 2006. 

Combining genetic and geospatial analyses to infer 

population extinction in mygalomorph spiders endemic to 

the Los Angeles region. Anim. Conserv. 9:145-157. 
Bray, Y., S. Devillard, E. Marboutin, B. Mauvy, and R. Peroux. 

2007. Natal dispersal of European hare in France. J. Zool. 

273:426-434. 

Brookfield, J. F. Y. 1996. A simple new method for estimating 

null allele frequency from heterozygote deficiency. Mol. 

Ecol. 5:453-455. 
Brooks, R. T. 2003. Abundance, distribution, trends, and 

ownership patterns of early successional forests in the 

northeastern United States. For. Ecol. Manage. 185:65-74. 
Brubaker, D. R. 2012. Monitoring occupancy and abundance 

of New England cottontails using noninvasive genetic tools. 

[Masters thesis], Univ. of New Hampshire, Durham, NH. p. 

129. 

Brubaker, D. R., A. I. Kovach, M. J. Ducey, W. J. Jakubas, and 
K. M. O'Brien. 2014. Factors influencing detection in 
occupancy surveys of an endangered lagomorph. Wildl. Soc. 
Bull, doi: 10.1002/wsb.416. 

Chambers, J. L., and D. Garant. 2010. Determinants of 

population genetic structure in eastern chipmunks (Tamias 
striatus): the role of landscape barriers and sex-biased 
dispersal. J. Hered. 101:413-422. 

Chapman, J. A. 1975. Syivilagus transitionalis. Mamm. Species 
55:1-4. 

Chapman, J. A., and G. G. Ceballos. 1990. The cottontails. Pp. 
95-110 in J. A. Chapman and J. E. C. Flux, eds. Rabbits, 
hares and pikas. I.U.C.N., Gland, Switzerland, p. 168. 

Chapman, J. A., and G. A. Feldhamer. 1981. Syivilagus 
aquaticus. Mamm. Species 151:1-4. 

Chapman, J. A., and G. R. Willner. 1981. Syivilagus paiustris. 
Mamm. Species 153:1-3. 

Chapuis, M. P., and A. Estoup. 2007. Microsatellite null alleles 
and estimation of population differentiation. Mol. Biol. 
Evol. 24:621-631. 

Chen, C, E. Durand, F. Forbes, and O. Francois. 2007. 
Bayesian clustering algorithms ascertaining spatial 
population structure: a new computer program and a 
comparison study. Mol. Ecol. Notes 7:747-756. 

Chevan, A., and M. Sutherland. 1991. Hierarchical 
partitioning. Am. Stat. 45:90-96. 

Clevenger, A. P., B. Chruszcz, and K. Gunson. 2001. Drainage 
culverts as habitat linkages and factors affecting passage by 
mammals. J. Appl. Ecol. 38:1340-1349. 

Coulon, A., G. Guillot, J.-F. Cosson, M. A. Angibault, S. 
Aulagnier, B. Cargnelutti, et al. 2006. Genetic structure is 
influenced by landscape features: empirical evidence from a 
roe deer population. Mol. Ecol. 15:1669-1679. 



1870 



© 2014 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



L. E. Fenderson ef al. 



Gene Flow in the New England Cottontail 



Crispo, E., J.-S. Moore, J. A. Lee-Yaw, S. M. Gray, and B. C. 
Haller. 2011. Broken barriers: human-induced changes to 
gene flow and introgression in animals. BioEssays 33:508- 
518. 

Crooks, K. R., and M. Sanjayan. 2006. Connectivity 
conservation. Cambridge Univ. Press, New York. 

Cullingham, C. L, C. J. Kyle, B. A. Pond, E. E. Rees, and B. N. 
White. 2009. Differential permeability of rivers to raccoon 
gene flow corresponds to rabies incidence in Ontario, 
Canada. Mol. Ecol. 18:43-53. 

Cunningham, K. M., M. F. Canino, I. B. Spies, and L. Hauser. 
2009. Genetic isolation by distance and localized fjord 
population structure in Pacific cod (Gadus macrocephalus): 
limited effective dispersal in the northeastern Pacific Ocean. 
Can. J. Fish. Aquat. Sci. 66:153-166. 

Cushman, S. A., and E. L. Landguth. 2010. Scale dependent 
inference in landscape genetics. Landscape Ecol. 25:967-979. 

Diffendorfer, J. E., and N. A. Slade. 2002. Long-distance 
movements in cotton rats (Sigmodon hispidus) and prairie 
voles (Microtus ochrogaster) in northeastern Kansas. Am. 
Midi. Nat. 148:309-319. 

Dixon, J. D., M. K. Oli, M. C. Wooten, T. H. Eason, J. W. 
McCown, and M. W. Cunningham. 2007. Genetic 
consequences of habitat fragmentation and loss: the case of 
the Florida black bear (Ursus americanus floridanus) . 
Conserv. Genet. 8:455-464. 

Dodd, C. K. Jr, W. J. Barichivich, and L. L. Smith. 2004. 
Effectiveness of a barrier wall and culverts in reducing 
wildlife mortality on a heavily traveled highway in Florida. 
Biol. Conserv. 118:619-631. 

Durand, E., C. Chen, and O. Francois. 2009. Tess version 2.3 
- Reference Manual [Documentation file, online]. Document 
last modified in August, 2009. Available at http://membres 
-timc.imag.fr/01ivier.Francois/manual.pdf. (accessed 10 
February 2014). 

Estes-Zumpf, W. A., and J. L. Rachlow. 2009. Natal dispersal 
by the pygmy rabbit (Brachylagus idahoensis) . J. Mammal. 
90:363-372. 

Estes-Zumpf, W. A., J. L. Rachlow, L. P. Waits, and K. I. 

Warheit. 2010. Dispersal, gene flow, and population genetic 

structure in the pygmy rabbit (Brachylagus idahoensis). 

J. Mammal. 91:208-219. 
Etherington, T. R. 2011. Python based GIS tools for landscape 

genetics: visualising genetic relatedness and measuring 

landscape connectivity. Methods Ecol. Evol. 2:52-55. 
Ewers, R. M., and R. K. Didham. 2006. Confounding factors 

in the detection of species responses to habitat 

fragmentation. Biol. Rev. 81:117-142. 
Fahrig, L. 2003. Effects of habitat fragmentation on 

biodiversity. Annu. Rev. Ecol. Syst. 34:487-515. 
Fahrig, L., and T. Rytwinski. 2009. Effects of roads on animal 

abundance: an empirical review and synthesis. Ecol. Soc. 

14:21. 



Fenderson, L. E. 2010. Landscape genetics of the New England 
cottontail: effects of habitat fragmentation population 
genetic structure and dispersal. [MS thesis], Univ. of New 
Hampshire, Durham, NH, p. 184. 

Fenderson, L. E., A. I. Kovach, J. A. Litvaitis, and M. K. 
Litvaitis. 2011. Population genetic structure and history of 
fragmented remnant populations of the New England 
cottontail (Sylvilagus transitionalis) . Conserv. Genet. 12: 
943-958. 

Forman, R. T. T., D. Sperling, J. A. Bissonette, A. P. 
Clevenger, C. D. Cutshall, V. H. Dale, et al. 2003. Road 
ecology, science and solutions. Island Press, Washington, 
DC. 

Frankham, R. 1995. Effective population size/adult population 
size ratios in wildlife: a review. Genet. Res. 66:95-107. 

Frankham, R. 2005. Genetics and extinction. Biol. Conserv. 
126:131-140. 

Frankham, R., J. D. Ballou, and D. A. Briscoe. 2002. 

Introduction to conservation genetics. Cambridge Univ. 

Press, New York. 
Franklin, I. 1980. Evolutionary change in small populations. 

Pp. 135-149 in M. E. Soule and B. A. Wilcox, eds. 

Conservation biology: an evolutionary ecological perspective. 

Sinauer Associates, Sunderland. 
Franklin, I. R., and R. Frankham. 1998. How large must 

populations be to retain evolutionary potential? Anim. 

Conserv. 1:69-70. 
Frantz, A. C, L. C. Pope, P. J. Carpenter, T. J. Roper, G. J. 

Wilson, R. J. Delahay, et al. 2003. Reliable microsatellite 

genotyping of the Eurasian badger (Meles meles) using faecal 

DNA. Mol. Ecol. 12:1649-1661. 
Frantz, A. C, L. C. Pope, T. R. Etherington, G. J. Wilson, and 

T. Burke. 2010. Using isolation-by-distance-based 

approaches to assess the barrier effect of linear landscape 

elements on badger (Meles meles) dispersal. Mol. Ecol. 

19:1663-1674. 

Fuller, S., and A. Tur. 2012. Conservation strategy for the New 
England cottontail (Sylvilagus transitionalis). Regional New 
England Cottontail Initiative, p. 148. Available at www. 
newenglandcottontail.org. (accessed 15 July 2013). 

Garza, J. C, and E. G. Williamson. 2001. Detection of 

reduction in population size using data from microsatellite 
loci. Mol. Ecol. 10:305-318. 

Gauffre, B., A. Estoup, V. Bretagnolle, and J. F. Cosson. 2008. 
Spatial genetic structure of a small rodent in a 
heterogeneous landscape. Mol. Ecol. 17:4619-4629. 

Gillis, E. A., and C. J. Krebs. 1999. Natal dispersal of snowshoe 
hares during a cyclic population increase. J. Mammal. 
80:933-939. 

Gonzalez, A., J. H. Lawton, F. S. Gilbert, T. M. Blackburn, and 
I. Evans-Freke. 1998. Metapopulation dynamics, abundance, 
and distribution in a microecosystem. Science 281:2045- 
2047. 



2014 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



1871 



Gene Flow in the New England Cottontail 



L. E. Fenderson ef al. 



Goslee, S. C, and D. L. Urban. 2007. The ecodist package for 
dissimilarity-based analysis of ecological data. J. Stat. Softw. 
22:1-19. 

Goudet, J. 1995. FSTAT (Version 1.2): a computer program to 

calculate F-statistics. J. Hered. 86:485^86. 
Graves, T. A., P. Beier, and J. A. Royle. 2013. Current 

approaches using genetic distances produce poor estimates 

of landscape resistance to interindividual dispersal. Mol. 

Ecol. 22:3888-3903. 
Guillot, G., and F. Rousset. 2013. Dismanteling the Mantel 

tests. Methods Ecol. Evol. 4: 336-344. 
Hardy, O. J., and X. Vekemans. 2002. SPAGeDi: a versatile 

computer program to analyse spatial genetic structure at the 

individual or population levels. Mol. Ecol. Notes 2:618-620. 
Holderegger, R., and M. Di Giulio. 2010. The genetic effects of 

roads: a review of empirical evidence. Basic Appl. Ecol. 

11:522-531. 

Jackson, S. N. 1973. Distribution of cottontail rabbits 

(Sylvilagus spp.) in northern New England. [M.S. thesis]. 

Univ. of Connecticut, Storrs, Connecticut. 
Jakobsson, M., and N. A. Rosenberg. 2007. CLUMPP: a cluster 

matching and permutation program for dealing with label 

switching and multimodality in analysis of population 

structure. Bioinformatics 23:1801-1806. 
Jamieson, I. G., and F. W. Allendorf. 2012. How does the 50/ 

500 rule apply to MVPs? Trends Ecol. Evol. 27:578-584. 
Jaquiery, J., T. Broquet, A. H. Hirzel, J. Yearsley, and N. 

Perrin. 2011. Inferring landscape effects on dispersal from 

genetic distances: how far can we go? Mol. Ecol. 20:692-705. 
Kalinowski, S. T. 2005. HP-RARE 1.0: a computer program for 

performing rarefaction on measures of allelic richness. Mol. 

Ecol. Notes 5:187-189. 
Keyghobadi, N. 2007. The genetic implications of habitat 

fragmentation for animals. Can. J. Zool. 85:1049-1064. 
Keyghobadi, N., J. Roland, and C. Strobeck. 1999. Influence of 

landscape on the population genetic structure of the alpine 

butterfly Parnassius smintheus (Papilionidae). Mol. Ecol. 

8:1481-1495. 

Kilpatrick, H., T. Goodie, and A. I. Kovach. 2013. Comparison 
of live-trapping and noninvasive genetic sampling to assess 
patch occupancy by New England cottontail rabbits. Wildl. 
Soc. Bull. 37:901-905. 

Kovach, A. I., M. K. Litvaitis, and J. A. Litvaitis. 2003. 
Evaluation of fecal mtDNA analysis as a method to 
determine the geographic distribution of a rare lagomorph. 
Wildl. Soc. Bull. 31:1061-1065. 

Laurence, S., M. J. Smith, and A. I. Schulte-Hostedde. 2013. 
Effects of structural connectivity on fine scale population 
genetic structure of muskrat, Ondatra zibethicus. Ecol. Evol. 
3:3524-3535. 

Legendre, P., and M. E. Fortin. 2010. Comparison of the 
Mantel test and alternative approaches for detecting 
complex multivariate relationships in the spatial analysis of 
genetic data. Mol. Ecol. Res. 10:831-844. 



1872 



Lichstein, J. W. 2007. Multiple regression on distance matrices: 
a multivariate spatial analysis tool. Plant Ecol. 188:117-131. 

Lindenmayer, D., R. J. Hobbs, R. Montague-Drake, J. 

Alexandra, A. Bennett, M. Burgman, et al. 2008. A checklist 
for ecological management of landscapes for conservation. 
Ecol. Lett. 11:78-91. 

Litvaitis, J. A. 2003. Are pre-Columbian conditions relevant 
baselines for managed forests in the northeastern United 
States? For. Ecol. Manage. 185:113-126. 

Litvaitis, J. A., and W. J. Jakubas. 2004. New England 
cottontail (Sylvilagus transitionalis) assessment. Maine 
Department of Inland Fisheries and Wildlife, Augusta, ME. 
p. 73. 

Litvaitis, M. K., and J. A. Litvaitis. 1996. Using mitochondrial 
DNA to inventoiy the distribution of remnant populations 
of New England cottontails. Wildl. Soc. Bull. 24:725-730. 

Litvaitis, J. A., and R. Villafuerte. 1996. Factors affecting the 
persistence of New England cottontail metapopulations: 
the role of habitat management. Wildl. Soc. Bull. 24:686- 
693. 

Litvaitis, J. A., B. Johnson, W. Jakubas, and K. Morris. 2003. 
Distribution and habitat features associated with remnant 
populations of New England cottontails in Maine. Can. J. 
Zool. 81:877-887. 

Litvaitis, J. A., J. P. Tash, M. K. Litvaitis, M. N. Marchand, A. 
I. Kovach, and R. Innes. 2006. A range-wide survey to 
determine the current distribution of New England 
cottontails. Wildl. Soc. Bull. 34:1190-1197. 

Litvaitis, J. A., M. S. Barbour, A. L. Brown, A. I. Kovach, J. D. 
Oehler, B. L. Probert, et al. 2008. Testing multiple 
hypotheses to identify causes of the decline of a lagomorph 
species: the New England cottontail as a case study. Pp. 
167-185 in P. Alves and K. Hacklander, eds. Biology of 
lagomorphs-evolution, ecology and conservation. Springer, 
New York. 

Lorimer, C. G., and A. S. White. 2003. Scale and frequency of 
natural disturbances in the northeastern US: implications 
for early successional forest habitats and regional age 
distributions. For. Ecol. Manage. 185:41-64. 

Lowe, W. H., and M. A. McPeek. 2012. Can natural selection 
maintain long-distance dispersal? Insight from a stream 
salamander system. Ecol. Evol. 26:11-24. 

Luikart, G., F. W. Allendorf, J. M. Cornuet, and W. B. 
Sherwin. 1998. Distortion of allele frequency distributions 
provides a test for recent population bottlenecks. J. Hered. 
89:238-247. 

Manly, B. F. J. 1986. Multivariate statistical methods. A 

primer. Chapman and Hall, London. 
Mantel, N. 1967. The detection of disease clustering and a 

generalized regression approach. Cancer Res. 27:209-220. 
Mata, C, I. Hervas, J. Herranz, F. Suarez, and J. E. Malo. 

2008. Are motorway wildlife passages worth building? 

Vertebrate use of road-crossing structures on a Spanish 

motorway. J. Environ. Manage. 88:407-415. 



© 2014 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



L. E. Fenderson ef al. 



Gene Flow in the New England Cottontail 



McKelvey, K. S., and M. K. Schwartz. 2005. DROPOUT: a 
program to identify problem loci and samples for 
noninvasive genetic samples in a cap-mark-recapture 
framework. Mol. Ecol. Notes 5:716-718. 

McRae, B. H., P. Beier, L. E. Dewald, L. Y. Huynh, and 
P. Keim. 2005. Habitat barriers limit gene flow and 
illuminate historical events in a wide-ranging carnivore, the 
American puma. Mol. Ecol. 14:1965-1977. 

MDIFW. 2007. "Recommended changes to Maine's list of 
endangered and threatened species." Maine Department of 
Inland Fisheries and Wildlife. Available at http://www. 
maine.gov/ifw/wildlife/species/pdfs/etlist_recommendations. 
pdf (accessed 5 April 2014). 

Newman, D., and D. Pilson. 1997. Increased probability of 
extinction due to decreased genetic effective population size: 
experimental populations of Clarkia pulchella. Evolution 
51:354-362. 

Newmark, W. D. 1987. A land-bridge island perspective on 

mammalian extinctions in western North American parks. 

Nature 325:430-432. 
Newmark, W. D. 1995. Extinction of mammal populations in 

western North American national parks. Conserv. Biol. 

9:512-526. 

NHFG 2008. "Endangered and threatened wildlife of New 
Hampshire" New Hampshire Fish and Game. Available 
from http://www.wildlife.state.nh.us/Wildlife/Nongame/ 
Nongame_PDFs/Endangered_Threatened_Wildlife_NH_ 
1108.pdf. (accessed 5 April 2014). 

Paetkau, D., R. Slade, M. Burden, and A. Estoup. 2004. 
Genetic assignment methods for the direct, real-time 
estimation of migration rate: a simulation-based exploration 
of accuracy and power. Mol. Ecol. 13:55-65. 

Perez-Espona, S., F. J. Perez-Barben'a, J. E. McLeod, D. Jiggins, 
I. J. Gordon, and J. M. Pemberton. 2008. Landscape features 
affect gene flow of Scottish Highland red deer (Cervus 
elaphus). Mol. Ecol. 17:981-996. 

Pimm, S. L., H. L. Jones, and J. Diamond. 1988. On the risk 
of extinction. Am. Nat. 132:757-785. 

Piry, S., G. Luikart, and J. M. Cornuet. 1999. BOTTLENECK: 
a computer program for detecting recent reductions in the 
effective population size using allele frequency data. 
J. Hered. 90:502-503. 

Piry, S., A. Alapetite, J. M. Cornuet, D. Paetkau, L. Baudoin, 
and A. Estoup. 2004. GENECLASS2: a software for genetic 
assignment and first-generation migrant detection. J. Hered. 
95:536-539. 

Pompanon, F., A. Bonin, E. Bellemain, and P. Taberlet. 2005. 

Genotyping errors: causes, consequences and solutions. Nat. 

Rev. Genet. 6:847-859. 
Pritchard, J. K., M. Stephens, and P. Donnelly. 2000. Inference 

of population structure using multilocus genotype data. 

Genetics 155:945-959. 
R Development Core Team 2006. R: a language and 

environment for statistical computing. R Development Core 



Team, Vienna, Austria. Available at http://www.R-project. 
org. (accessed 10 December 2010). 
Rannala, B., and J. L. Mountain. 1997. Detecting immigration 
by using multilocus genotypes. Proc. Natl Acad. Sci. USA 
94:9197-9201. 

Raymond, M., and F. Rousset. 1995. GENEPOP (version 1.2): 

population genetics software for exact tests and 

ecumenicism. J. Hered. 86:248-249. 
Reed, D. H., J. J. O'Grady, B. W. Brook, J. D. Ballou, and 

R. Frankham. 2003. Estimates of minimum viable 

population sizes for vertebrates and factors influencing those 

estimates. Biol. Conserv. 113:23-34. 
Reed, D. H., A. C. Nichols, and G. E. Stratton. 2007. Genetic 

quality of individuals impacts population dynamics. Anim. 

Conserv. 10:275-283. 
Rosenberg, N. A. 2004. Distruct: a program for the graphical 

display of population structure. Mol. Ecol. Notes 4:137- 

138. 

Rothermel, B. B., and R. D. Semlitsch. 2002. An experimental 
investigation of landscape resistance of forest versus 
old-field habitats to emigrating juvenile amphibians. 
Conserv. Biol. 16:1324-1332. 

Rousset, F. 2000. Genetic differentiation between individuals. 
J. Evol. Biol. 13:58-62. 

Sauer, J. R., J. E. Hines, J. E. Fallon, K. L. Pardieck, D. J. 
Ziolkowski Jr, and W. A. Link. 2011. The North American 
breeding bird survey, results and analysis 1966-2009. 
Version 3.23.2011. USGS Patuxent Wildlife Research Center, 
Laurel, MD. 

Schwartz, M. K., and K. S. McKelvey. 2009. Why sampling 
scheme matters: the effect of sampling scheme on landscape 
genetic results. Conserv. Genet. 10:441-452. 

Schwartz, M. K., L. S. Mills, Y. Ortega, L. F. Ruggiero, and F. 
W. Allendorf. 2003. Landscape location affects genetic 
variation of Canada lynx. Mol. Ecol. 12:1807-1816. 

Segelbacher, G., J. Hoglund, and I. Storch. 2003. From 

connectivity to isolation: genetic consequences of population 
fragmentation in capercaillie across Europe. Mol. Ecol. 
12:1773-1780. 

Segelbacher, G., S. A. Cushman, B. K. Epperson, M. Fortin, O. 

Francois, O. J. Hardy, et al. 2010. Applications of landscape 

genetics in conservation biology: concepts and challenges. 

Conserv. Genet. 11:375-385. 
Storfer, A., M. A. Murphy, J. S. Evans, C. S. Goldberg, S. 

Robinson, S. F. Spear, et al. 2007. Putting the 'landscape' in 

landscape genetics. Heredity 98:128-142. 
Taberlet, P., S. Griffin, B. Goossens, S. Questiau, V. Manceau, 

N. Escaravage, et al. 1996. Reliable genotyping of samples 

with very low DNA quantities using PCR. Nucleic Acids 

Res. 24:3189-3194. 
Tallmon, D. A., A. Koyuk, G. Luikart, and M. A. Beaumont. 

2008. ONeSAMP: a program to estimate effective 

population size using approximate Bayesian computation. 

Mol. Ecol. Res. 8:299-301. 



2014 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



1873 



Gene Flow in the New England Cottontail 



L. E. Fenderson ef al. 



Tash, J. P., and J. A. Litvaitis. 2007. Characteristics of 

occupied habitats and identification of sites for restoration 
and translocation of New England cottontail populations. 
Biol. Conserv. 137:584-598. 

Templeton, A. R., K. Shaw, E. Routman, and S. K. Davis. 
1990. The genetic consequences of habitat fragmentation. 
Ann. Mo. Bot. Card. 77:13-27. 

Templeton, A. R., R. J. Robertson, J. Brisson, and J. Strasburg. 
2001. Disrupting evolutionary processes: The effect of 
habitat fragmentation on collared lizards in the Missouri 
Ozarks. Proc. Natl Acad. Sci. 98:5426-5432. 

Tischendorf, L., and L. Fahrig. 2000. On the usage and 
measurement of landscape connectivity. Oikos 90:7-19. 

USFWS. 2006. Endangered and threatened wildlife and plants; 
review of native species that are candidates or proposed for 
listing as endangered or threatened; annual notice of 
findings on resubmitted petitions; annual description of 
progress on listing actions. United States Fish and Wildlife 
Service. Fed. Reg. 71:53755-53835. 

USFWS. 2012. Endangered and threatened wildlife and plants; 
review of native species that are candidates or proposed for 
listing as endangered or threatened; annual notice of 
findings on resubmitted petitions; annual description of 
progress on listing actions. United States Fish and Wildlife 
Service. Fed. Reg. 77:69994-70060. 

Van Oosterhout, C, W. F. Hutchinson, D. P. M. Wills, and P. 
Shipley. 2004. Microchecker: software for identifying and 
correcting genotyping errors in microsatellite data. Mol. 
Ecol. Notes 4:535-538. 

Verhoeven, K. J. F., K. L. Simonsen, and L. M. Mclntyre. 
2005. Implementing false discovery rate control: increasing 
your power. Oikos 108:643-647. 

Walsh, C, and R. MacNally. 2008. Package 'hier.part': 
hierarchical partitioning. R package version 1.0-3. R 
Foundation for Statistical Computing, Vienna, Austria. 

Waples, R. S. 2006. A bias correction for estimates of effective 
population size based on linkage disequilibrium at unlinked 
gene loci. Conserv. Genet. 7:167-184. 

Willi, Y., J. van Buskirk, B. Schmid, and M. Fischer. 2007. 
Genetic isolation of fragmented populations is 
exacerbated by drift and selection. J. Evol. Biol. 20: 
534-542. 

Williamson-Nateson, E. G. 2005. Comparison of methods for 
detecting bottlenecks from microsatellite loci. Conserv. 
Genet. 6:551-562. 

Zalewski, A., S. B. Piertney, H. Zalewska, and X. Lambin. 
2009. Landscape barriers reduce gene flow in an invasive 
carnivore: geographical and local genetic structure of 
American mink in Scotland. Mol. Ecol. 18:1601-1616. 

Zartman, C. E., S. F. McDaniel, and A. J. Shaw. 2006. 
Experimental habitat fragmentation increases linkage 
disequilibrium but does not affect genetic diversity or 
population structure in the Amazonian liverwort Radula 
flaccid. Mol. Ecol. 15:2305-2315. 



1874 



Appendix 1: GIS Data Sources and 
General Processing Methods 

Roads Model 
Data source 

• U.S. Geological Survey (USGS), Maine Office of Geo- 
graphic Information Systems (MEGIS) (ed.) trans, 
1989: Augusta, Maine. Available online at http://megis. 
maine.gov/catalog. Accessed 3 May 2010. 

We obtained the "trans" feature dataset from Maine GIS 
(which extends into the small area of New Hampshire 
included in our dataset) and clipped it to a rectangle outlining 
the extent of our sampling locations. This file includes roads, 
trails, pipelines, railroads, and powerline utility corridors 
("Otrans") to be used in the development of the facilitator 
models. To create a vector file of the roads, we first selected the 
class 1-6 roads and buffered the polyline roads by class, such 
that busier roads had a wider buffer, which would translate 
into wider polygons for interstates versus trails, for example. 
Buffer sizes were as follows: roads classified as "interstate" were 
buffered to 30 m, primary roads 20 m, secondary roads 10 m, 
improved roads 5 m, unimproved roads 2 m, and trails 1 m. 

Surficial Water Model 
Data source 

• U.S. Geological Survey (USGS) and U.S. Environmen- 
tal Protection Agency (USEPA). National Hydrography 
Dataset (NHD) Medium Resolution, 2005; updated in 
2006 with SOSC (Augmenting NHDPlus Strahler order 
values using Strahler calculator): Corvallis, Oregon. 
Available online at http://nhd.usgs.gov. Accessed 25 
March 2011. 

To create a model of surficial water, we first merged the 
NHDwaterbody (lakes, ponds, etc.) and NHDArea (rivers, 
inlets, etc.) files. This output was then joined with the 
NHDflowline shapefile, based on spatial location, such that 
the maximum Strahler stream order class was assigned to 
each water feature. Water features without stream order 
information were given a stream order of 1. 

Facilitator Model 

Data sources 

• U.S. Geological Survey (USGS), Maine Office of Geo- 
graphic Information Systems (MEGIS) (ed.) trans, 
1989: Augusta, Maine. Available online at http://megis. 
maine.gov/catalog. Accessed 3 May 2010. 
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• Roads. Created by Lindsey Fenderson, as a modifica- 
tion of trans, using ArcGIS 9.3, as described above. 
(May 2010). 

Road edges were established by buffering class 1-2 roads 
from the above-created "Roads" file by 30 meters and buf- 
fering class 3-6 roads by 10 m (as verges maintained on 
highways are typically wider than those found on reduced- 
volume roads), then erasing the roads themselves. Powerline 
and railroad features were obtained from the "Otrans" of 



the clipped "trans" dataset. These polylines were buffered by 
30 m and combined with the "RdEdges" shapefile for the 
final "Facilitators" dataset. 

All vector models (roads, surficial water, and facilitators) 
were then converted to rasters with a cell size of 10 for cost- 
distance analyses. Raster files were reclassed according to 
values given in Table 2 to develop each cost model per fea- 
ture type (binary, linear, exponential, and logarithmic). 
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